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We present a systematic quasi-mean field model of the Ostwald ripening process in two dimen- 
sions. Our approach yields a set of dynamic equations for the temporal evolution of the minority 
phase droplets' radii. The equations contain only pairwise interactions between the droplets; these 
interactions are evaluated in a mean- field type manner. We proceed to solve numerically the dy- 
namic equations for systems of tens of thousands of interacting droplets. The numerical results are 
compared with the experimental data obtained by Krichevsky and Stavans for the relatively large 
volume fraction tp — 0.13. We found good agreement with experiment even for various correlation 
functions. 



I. INTRODUCTION 

When a system is quenched into a two-phase coexistence region, its homogeneous initial state no longer corresponds 
to thermodynamic equilibrium!!!. If there is a conserved quantity whose density is different in the two coexisting phases 
(such as the total volume or total amount of impurities), the final equilibrium state consists of two macroscopic 
domains, separated by a single phase boundary. The manner in which a system evolves from its homogeneous 
initial state to the final equilibrium state of two-phase coexistence has been the subject of numerous theoretical and 
experimental investigations. 

For small initial super-saturation the system's evolution towards two-phase equilibrium starts by nucleation and 
growth of droplets (or small crystallites) of the minority phase. In the late stage of evolution to the equilibrium state 
no new droplets are formed and the amount of material in each of the phases remains fixed. Evolution proceeds by 
means of dissolution of small droplets and growth of the large ones, giving rise to reduction of the total (surface) energy 
of the system. The exchange of material is driven by diffusion; the concentration is higher near surfaces with high 
curvature, and hence the diffusive flux is directed from the small droplets towards the larger ones. This coarsening 
processB, called Ostwald ripening, in the course of which the number of droplets decreases, while the average size of 
the remaining ones increases, is the subject of the present paper. 

In what follows we describe a formalism that leads to an efficient numerical algorithm for Ostwald ripening in 
two dimensions. Our work was motivated by, and our results are compared with recent experimental work on a 
two-dimensional film of liquid and crystalline succinonitrile in coexistenceu. Therefore we will refer to the minority 
phase as "solid", and to the matrix in which the droplets are embedded as "liquid". 

Ostwald ripening belongs to a family of non-equilibrium phenomena that exhibit evolution to a scaling state. By a 
scaling state we mean that there is a single length scale in the system, which grows with time. Forjpstwald ripening 
the obvious length scale is the average droplet size (R) ; its growth with time follows the celebratedcl Lifshitz-Slyozov 
law, (R) ~ t 1 / 3 . When rescaled by this changing length, all statistical characteristics of the system (such as droplet 
size distribution, spatial correlations etc) are time-independent. That is, by comparing two photographs of the system, 
taken at different (well separated) times but properly rescaled, one cannot tell from any statistical measurement which 
photograph was taken earlier. 

A parameter of central importance on which various characteristics of the system (such as the droplet size distri- 
bution) depend is the relative volume fraction of the minority phase, p. This parameter determines the ratio of the 
size of the droplets and the distance between them. Therefore, p controls the extent to which droplets interact with 
each other; for very small volume fractions the interaction is weak, and any particular droplet "feels" the effect of 
the othees only through an effective medium. This observation served as the basis of the theoretical, mean - field 
modelsoLl (see ref.B for a review). As ip increases correlations between neighboring droplets become more important, 
until mean field based approximations lose their, validity. The effects of correlations have been taken into account 
analytically to first order in a small parameteitHl3. In three dimensions the small parameter is ^JTp (see Mardem) and 
in two - l/log<^ (see Zheng and GuntontlJ). The obtained expressions are so complicated even in first order, that it 
is hard to see how one can proceed to higher orders. Rather, a combined analytical and numerical approach seems to 
be suitable. 

The basic phenopaenological description of Ostwald ripening is provided in the framework of a Cahn-Hilliard 
equatiorilJ (see ref.liHl for a review) for the order parameter. Rogers and Desait 3 ! performed such a "first princi- 
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pie" computation, reaching the scaling state with about 500 droplets. This was achieved for small values of <p (less 
than 0.1), when all the droplets were almost circular. Such a number of droplets suffices for studying the droplet size 
distribution, but is too small to gain insight into the spatial and temporal correlations in the scaling state. Moreover, 
for larger fractions the final stage of the coarsening process was not achieved in the calculation and scaling was not 
observed. MasbaumEJ used this method recently for a simulation, performed on a parallel computer, starting with 
about 3000 droplets. He obtained rather good agreement for various correlation functions with the experiments of 
Krichevsky and StavansB. 

It is clear, however, that this method can hardly be applied for studying much larger systems, which are needed in 
order to make statistically meaningful measurements in the late-stage scaling regime. Therefore one turns to simplified 
descriptions of the system, hoping that the main features of the full Cahn-Hilliard theory are preserved. For the sake 
of completeness, we briefly present below a sequence of approximations that lead from the full Cahn-Hilliard theory 
to our approach. 

First, one goes to a coarse-grained representation of the order parameter field, which retains the boundaries between 
the crystalline droplets and the surrounding liquid as one of the dynamic variables of the problem. The second variable 
is the diffusing concentration field (of the impurities or the liquid itself), c(r). Simplification of the problem is attained 
by separation of time scales. The process which occurs on fastest scales is the equilibration of the concentration near 
the (moving) boundaries. Assuming that this occurs instantaneously allows us to write the Gibbs-Thomson condition 
for the concentration at a point on the boundary of a droplet, where the (local) radius of curvature is R(t); 

OL 

c\droplet(t) = C eq (R(t)) — Coo + 77777' (!) 

H[t) 

Here Coo is the equilibrium concentration in a liquid above a planar liquid-solid interface. The problem becomes now 
one of solving the diffusion equation 

with boundary conditions given by Eq.(^) on moving boundaries. 

The rate of change of the boundaries is determined by Jl = n ■ J, the mass flux normal to the droplet's surface, 
which, in turn, is proportional to the gradient of the concentration field at the boundary; 

f = -J± J = - Vc k (3) 

Equations Eq.^l) constitute a closed dynamical problem. The next simplification is a quasi-static approximation 
to the problem. For fixed boundary conditions the diffusing field would reach a steady state in a characteristic diffusion 
time trj ~ R 2 , where R is a typical length scale (distance between neighboring droplets). If this time is much shorter 
than the growth time to, i.e. the time it takes a typical droplet's radius to change appreciably, the concentration field 
reaches a steady state before the radii of the droplets had a chance to change and one can use the stationary diffusion 
equation instead of Eq.(||). This is precisely the case for late stages of the growth process, where to ~ R 3 , so that 
indeed to -C tc- The Lifshitz-Slyozov scaling law, tn ~ R 3 , can be established by dimensional analysis of Eqs.(Q) and 
(|J). Hence at late stages the solution c(r, t) of Eq.(pJ) can be approximated by the solution of the stationary diffusion 
problem 

V 2 c = (4) 

with the previous boundary conditions Eq.(Q); but now we look for a solution of Eq.(Q), from which the rate of change 
of the boundaries is obtained using Eq.(H), giving rise to new boundaries at the next time step, and so on. It should 
be noted that in this approximation the total area of the droplets is conserved exactly. 

It is possible to reformulate the stationary diffusion problem Eqs.(|l|) - (|[) as an integral equation which auto- 
matically accounts for the boundary conditions. Expanding the shapes of the droplets using a set of orthogonal 
polynomials one can-Deduce the problem to an implicit system of ordinary differential equations in terms of the 
expansion coefficicntsEjllZl. 

Further simplification can be achieved by neglecting the deviation of the droplet shapes from circular. We expect 
this approximation to work for low volume fractions ip, when the distances between the droplets are much larger than 
the droplets sizes and redistribution of material in a single droplet is much faster than the exchange between the well 
separated droplets. However, experimentsd show that even for fractions as large as <p = 0.4 the droplets are more 
or less circular. Therefore, we consider Laplace's equation Eq.(^) with the boundary conditions Eq.(|l|) at perfectly 
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circular domains of radii Ri, positioned at points r^. This is a problem in electrostatics - of calculating the potential 
in the presence of conducting cylinders. One approximates the solution of this problem by 

c(r) =$> log(|r - n\/Ro) + E V-*!^ ' (5) 

where i?o is an arbitrary length. The "charges" g< and the "dipoles" Pi should be chosen so that the boundary 
conditions ([!]) are satisfied (approximately) on the surface of each droplet. Let us denote by Xij = \fl — rj\ the 
distance between the centers of droplets i and j, and by c|_R i the correct boundary value of the concentration, as given 
by Eq.([j]). We now approximate the concentration field at the the boundaries by its expansion to first order in the 
small parameter R^/Xif. 

a - (v ■ R ) 

+ E ft m^/ho) + E + E ^ 

denotes the radius- vector of points on the surface of droplet i; the last term appears as the expansion of qj log \fi + 
Ri — rj\. Equation (||) contains two parts; one that depends on the direction of Ri and one that does not. This results 
in two sets of linear equations for the charges and the dipoles: 

Coc + £ = ft log Ri/Ro + E ft lGg(*« /i?0)+E ^y^^ > ( 7 ) 



■ Jfc) x (X tJ ■ Ri) 



^ " \X. 

As we will see, the last term in Eq.(Q) is 0((R/ X) 2 ) - i.e. smaller than the other terms and can be omitted. Thus 
the system of equations 

coo + = ft log R4/R0 + E ft lo s( x y / R o) (9) 
determines the charges. Since Eq.(||) should be satisfied for any R4, the dipoles are determined by 

| = -Eftr|%- do) 



The physical meaning of the Eq.(10) is that the dipolc associated with droplet i is simply related to the electric field 
induced at its center by all the other charges. Once the charges are obtained from solving Eq. (^), their substitution 
in Eq.(|l(i|) determines the dipoles. 

The normal component of the flux at the boundary of droplet i is given by 

V„c(r J + i? 4 ) = f -2^1^, (11) 

where Eq. (g) and (0) were used. The normal flux has two contributions; an isotropic part, giving rise to a rate of 
change of the radius, given by 

t - 1 <"> 



and an anisotropic part, due to the dipole term in Eq.(|llj). The contribution of the dipole part to the total flux 
vanishes; it induces deposition of material on one side of the droplet and evaporation from the opposite side. We 
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approximate the effect of the dipole flux by shifting the positions of the circular droplets (see below). We also show 
(in Appendix A) that in the low concentration (ip « 1) limit the effect of the dipoles is negligible; we first discuss 
this limit, working with charges only, and then include the dipoles. 

Eq^l2|) implies that the rate of change of the area of droplet i is proportional to its charge qi. Using qi from 
Eq.(|l2) in Eq.([)|) eliminates the concentration field c{r) from the problem and, after proper rescaling of variables, the 
following set of equations for the temporal evolution of the radii results: 

(l + Ri/Rc)=^L itj R j , (13) 

3 

where R c = a/coo is a capillary length and the matrix Lij defined in 2-d as followsEl: 

T -RR IMW if f = J ( , 

L, id - muj i ^(X^/Ho) otherwise. 1 > 

Clearly, solving the system Eq.(|p|)-(|l4|) is a much simpler computational task than solving the Laplace problem 
Eq.(^) with moving boundary conditions on the surfaces of the droplets. Note, however, that the matrix elements 
Lij grow with the distance between the droplets. Because of this, finite size effects become crucial: all droplets in 
the bulk feel the boundary. The simplest way to avoid this problem is to consider the system with periodic boundary 
conditions; the interaction o£a pair of droplets contains an infinite sum of logarithms, corresponding to all the images 
of these droplets. Yao et alia have summed this series using Ewald summation techniques. 

The problem is further complicated by the need to invert the N x N matrix Li j at each time step in order to 
get from Eq([l3"|) explicit expressions for R4. Thus, solving the system (p^-14]) takes N 3 operations per time step; 
therefore each run costs N 4 operations. This necessitates huge CPU times(the simulations of Yao et alE3 took about 
1000 CPU hours on an IBM 3090 computer for a system of about 3000 droplets). 

BeenakkerEj has solved an analogous problem in 3-d by truncating the matrix Li.j (defined in 3-d in a way analogous 
to Eq.(p^)) . He took into account only the interactions between droplets whose separation does not exceed a threshold, 
reducing this way the numher-of droplets with which a given one interacts to about 20. The physical motivation for 
such truncation is screening^Bfc}\; droplets whose separation exceeds the screening length do not affect each other. It 
should be noted that with such a truncation the total volume of the droplets is not conserved, and Beenakker had to 



adjust R c in Eq.(|14[)) at, each time step in order to restore volume conservation. 

Akaiwa and MeironllZl used the analogous truncation procedure in 2-d. Although the interactions between the 
charges are expected to be well screened even in 2-d, formal truncation of the matrix seems problematic in this 
case. Since the matrix elements grow with distance between the droplets, the elements of the inverse matrix as 
functions of the cutoff should contain large fast oscillating components (unlike in 3-d, where their dependence on the 
cutoff is smooth). Nevertheless, their results (in particular, the correlation functions) compare well with experiment. 
Apparently, this success is due to the fact that these oscillations are effectively averaged out during the run. 

Our goal was to find an analytic way of approximating L , in a manner that reflects screening as the physical 
basis of the approximation scheme. Once L has been inverted, we can integrate the equations 

ii, »-/ M I Hj R, I (15) 

3 

numerically. The approximation introduced in this work can be summarized by the expression 

o_y f -i_ 1 y Mx^/U ( 1 n f , 

1 2-< « RiK (Ri/Csc) K {R 3 /( SC ) \Ri RiJ' [ ' 

3 3\^r % ) 

where Kq is the zeroth order modified Bessel function of the second kind (also called MacDonald function). Since 
Kq(x) ~ exp (— x) for large x, the parameter £ sc has the meaning of a screening length. In the mean field limit (i.e. 
for very small area fraction ip) it is determined by the equation: 

c2=2 ™(5wk))- (17) 

However, for larger fractions, where the effects of correlations become important, £ sc is determined by a more com- 
plicated condition, as discussed in Appendix B. In practice, in our simulations for tp = 0.13 we used Csc — 2.13R, as 
determined by Table I (see Appendix B for details). 
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In spite of various approximations that were made to derive Eqs. (]T6|) - (17), we believe that they do contain the 
most important aspects of the dynamics present in Eq.(|l3|). One can see explicitly the effects of screening, with the 
screening length appearing naturally in the derivation presented below. Total area is conserved explicitly. In the 
limit of very small minority phase area fractions, where (i?)/£ sc — > 0, our description is consistent with Marqusee's 
mean- field theoryu, which leads to the following dynamical equation for the droplet 's radius: 

= WO (i - (18) 



dt v ' ^' \R R 

where k(x) = xK\(x) / Kq(x), while ( is defined by the equation: (~ 2 = 2wn (k(R/Q). Since xK±(x) — > 1 for x — > 0, 
we find k(x) — ► 1/Kq(x) and our definition of the screening length coincides with his. 

We used the formalism presented above to integrate the evolution of large assemblies of droplets and found that 
at relatively large values of the area fraction ip the droplets' motion, induced by the so far neglected dipoles becomes 
important and must be taken into account. 

We evaluate now the contribution of the dipoles to the shift of the droplets' centers of mass, 5fi, defined by 

M5r = J RSm. (19) 

Here the integration is over the boundary of a droplet; R is a radius vector on its boundary, M — pirR 2 is the total 
mass of the droplet (p is the density; it will drop out of the final result) and 5m is the additional mass adsorbed (or 
lost) locally, due to the shift of the boundary during the interval dt. The anisotropic component of the local velocity 
of the droplet 's boundary is given by 

2p cos 6 

V = - Jn = R? ' 

where 9 is the polar angle between R and p. The mass added at R is 

5m = p(vdt)(R dO) = -(2p/R)pdt cos 6 dJS, 

so that (see also Appendix A o%) 

df 2P f Rcos0 d6=-X [cos 2 0d8 2p 



dt irR 3 J ttR 2 J R 2 ' 

Finally, 



Note, that both sides of Eq.(p(i|) have the same dimensionality, since the charges qi are measured in area per time (see 
Eq.©). 

The procedure for solving Eq.(p^[)-(p^) for the dynamics of the droplets' radii Ri, together with Eq.(p0[) for their 
positions fj is as follows. For given and Ri we invert L, obtain dRi/dt (or qi as given by Eq([l2|)) and integrate one 
time step. Next, substituting q± into Eq.(|20|) we obtain new values of r*j and the procedure is repeated. 



The sum on the r.h.s. of Eq.(20) appears to be problematic: assuming that the droplets' charges are uncorrelated 
(while the total charge is zero) one can see that the mean square of this sum diverges logarithmically with the size 
of the system. However, as we show in Appendix A, the correlations between the charges, provided by screening, 
ensure the convergence of this sum, so that it can be evaluated using a reasonable fraction of its terms. Actually, a 
heuristic formula, simpler than that of Eq. (|20|) , can be used to evaluate the droplets' shifts in practical simulations. 
As explained in Appendix A, Eq.(pp[) can be replaced by 

^ = oS"^(l 1 ^ (21] 

dt ^ hki \Ri Ri) \XiA 2 ' 



where Ki j = K (Xij /(^ sc ) and fe, = Ko(Ri/£ sc ). Although the approximation leading to this formula is not clearly 
established, our simulations show that it works as well as the more rigorous Eq.(G(jh. At the same time it is much 
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more economic because only droplets that lie within the screening length (from the droplet whose shift is evaluated) 
contribute to the sum. 

In the sequence of approaches to the Ostwald ripening problem in two dimensions, ranging from the Cahn-Hilliard 
equation to Marqusee's mean field theory, our model, which includes only the pairwise interaction between the droplets, 
can be viewed as the minimal extension of the mean field approach. This concerns our Eq. ( pl| ) (in comparison with 



the Eq.(|20[)) as well. We believe that our approach constitutes the simplest step that can be taken beyond simple 
mean field. 

Our program is as follows. In the next section we present a few universal properties of the matrix L. This is 
followed by Section III, where we derive the main result, Eq^lq), using a mean-field typ e approximation to L^ 1 . 
Numerical solutions of the resulting dynamics, based on Eq.(|l6), and on both ( |20| ) or (plf), are presented in Section 
IV and compared with experiments by Stavans and Krichevsky. Section V presents a short summary of our approach 
and results. 

II. SOME SUM RULES 

We present now some important properties of the matrix L , that lead to useful "sum rules". The first claim is 
that 

^ /.,.;//, y,l.„:l<, •<• (22) 

i 3 

the sum approaching zexq as the system size increases. This relationship is analogous to one presented for the 3-d case 
by Beenakker and Ros£3. We present here arguments (that seem to us simpler than those of ref.t3) for the validity 
of Eq. (p2]) in 2-d. First consider the quantity 

*3 = E = E ^ i°g(^ARo). (23) 

i i 

The sum is clearly dominated by regions that are far from fj ; in the absence of long-range correlations the composition 
of such regions does not depend on the identity of droplet j. On the other hand, the contribution of those droplets 
i that lie near droplet j will, in general, depend on Rj. By "near" we mean the region for which correlations are 
important. The size of this region is, however, small (on the order of the screening area - see below), the values of 
Xij in this region are small, and the contribution from it is negligible compared to those from the far-away parts of 
the plane. That is, for increasing system size a 3 — ► oo, whereas the contribution from the correlation region remains 
finite. Hence we can write 

aj = Lij j Rj ~ a — > oo. (24) 

i 

Next, consider the following sum: 

i E^- EE L u L i* = E * -."j E L ^/ R i = E / -,; / ^ ( 25 ) 



using now Eq.([24|) we get 

"E^X- - 1 ( 26 ) 

3 

so that 

EAnM-^^O (27) 

3 

and hence statement Eq.(p2|) holds in the limit of large system size. An immediate consequence of Eq.(p7|) is obtained 
by using it in Eq.(|l5|), yielding 

E 7 -'.;- ( 28 ) 



G 



Thus we see that the length R c drops out from the description of the dynamics of the system. 
Multiplying by Ri and summing over i, and using again Eq.(p7j), we establish area conservation 



e^=e(e^) 



(29) 



as a consequence of our sum rule. Another important consequence is the following observation: , the elements of 
the inverse matrix, are independent of the parameter Rg. To see this, take the derivative of (L~ 1 L)i t k with respect to 
logi?o: clearly (L^ 1 )' L + L~ l L' = and hence 

(^i,m)' — ~ E^/^'*^™' (30) 

jk 

But from Eq.(|T4|) we immediately get that (Lj^)' — —RjRk, which when used in Eq.(p0[) gives 

( L iA)' = E E^.^'^./^ ( 31 ) 

3 k 

and when Eq.(^) is used here, we get 

(£,7™)' = ° ( 32 ) 

and hence LZ 1 does not depend on Rq; and in view of Eq.(|28|) the system's dynamics is, therefore, also independent 
ofi? - 

III. MEAN FIELD APPROXIMATION TO L _1 



Let us decompose the matrix L of Eq.(14) to its diagonal and off-diagonal parts, 

L = L -L U (33) 

where (L )ij — SijRf log(i?j/i? ) and = —RiRj\og(Xij/R ). Note, that only the off-diagonal part contains 

interactions between different droplets. One could consider an expansion of L _1 in powers of the interaction: 

oo 

L- 1 = (1 - L„ %) _1 £o 1 = E( Z o ^O" V (34) 



n=0 



with L Q 1 given by 



(Lo ^- 6 ^ Rnog{Rt/Ro y ( 35 ) 



Clearly, for matrix elements that correspond to well separated droplets ( i.e. when 3> Ro), terms of order n + 1 
will be larger than those of order n: that is, the scries diverges formally and in order to obtain meaningful results one 
should perform some kind of partial summation to all orders. To do this we introduce the T — matrix, defined by 

L- 1 =L 1 +L- 1 TL-\ (36) 
To ensure that Y2j -Li.jLJ^. — 8i t k, T has to satisfy the equation 

f = L 1 +L 1 L Q 1 f. (37) 
For the sake of convenience we introduce a new matrix 0, defined by 

Tij = RiRjftij. (38) 
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Once the matrix tfiij has been found, it is straightforward to write down T and, using Eq.(^), the inverse matrix 
L- 1 . Rewriting Eq.@ in terms of 4> we get 

For the diagonal elements this takes the form 

<t>k,k = - wp 7p log(Xkj/Ro)4>j,k, (40) 

whereas for the off-diagonal elements we obtain from Eq.(p9[) 

Yl 1 u~Td~ log(^ /i?o)^,fe + 4>i,k = -Ik logX, ik /Ro (41) 

where 

log(i?fe/i?o) 

It should be noted that equations (EclS) are exact. In principle for any configuration of droplets, the set of equations 
( |39l ) are to be solved for the matrix elements </>i,fe. We approximate the solution of this problem in a mean-field spirit. 
The central premise of this mean field approach is that the off-diagonal matrix elements 4>i,j depend only on the 
distance Xij, i.e. 

<^=^(H-r-|). (43) 

That is, we are interested only in the pairwise interaction between the droplets (neglecting any possible dependence 
of cf>i j on other droplets). The analogous approximation for the diagonal elements is that they are all equal, i.e. <j)™l 
is independent of k; 

<f>k,k « Cfc = (44) 
To obtain manageable mean-field equations the following simplifying assumptions are made: 

1. In sums such as ( (40f]4l| ) replace l/\og(Rj/Ro) by its average valueEl. 

2. Approximate these sums by integrals; 

3. Set in Eq. ( [fl| ) jk = J for all k (i.e. impose independence of k). 

We will discuss and justify these steps below, but before doing that, we investigate the resulting approximation to 
Eq.@, given by 

- ^Co 2 i log(X 1J /i?o)0(X J , fe )d 2 r J + 4>{X ifk ) = -jlogX itk /R , (45) 



where 



c " = 2 ™ {^km) ra 2 ™iogww) (46) 

and the angular brackets denote averaging over the distribution of droplet sizes (that is, the average value obtained 
in the particular droplet configuration in which L _1 is evaluated). 

The integral equation ( jl5| ) can be solved easily by operating with yf on the left and the right hand sides. Using 
the identity V 2 l°g r — 27T(5(r) we obtain the following differential equation for 4>(x): 

-C,rV + V 2 = -2^(rO- (47) 

Its solution is the MacDonald function, 
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(48) 



-log(r/2Co) -C r«Co 
K (r/Co)^{ r— j (49) 

i^- rKo r»Co 

where C « 0.5772 is Euler's constant. 

Let us discuss briefly this result. First of all note that self consistency of the approximation imposes positivity of 
Co, i-e. Rq ^S> (R). This is indeed the case, as will be discussed below. Next note that the divergence of <j){r) at short 
distances is an artifact of the approximations we made by replacing the exact discrete equation ( |4l| ) by the continuous 
Eq. ([f5"|) . This divergence has no physical consequences, since the diagonal elements of the discrete matrix (j) are not 
given by taking r — > in the solution of Eq.([47j); rather, they are determined by Eq. (f40|) . 

One should note also that if j k were to depend on k the result (Q) would have become </>i ;k — 7fc-ft"oPQ,fc/Co), 
implying that (f>i tk ^ <p k ,i, whereas the matrices T and (f> must be symmetric (see Eqs. ( |36| ) and ( ^8|) and remember 
that L is symmetric). 

The diagonal elements are expressed in terms of the off-diagonal ones in Eq. ( [40|); using there the approximations 
listed above and the mean-field expression (^) for the off-diagonal elements, Eq.([40]) becomes 



h.k « ^Co" / ^{Xkj/Ro)K {X j , k /<j a )d>r j = 0o (50) 

Z7T 



as anticipated in Eq.(MJ). Using this in Eq.(M3) yields 



7fe 



log(i? fc /i? ) 



and in order to get a consistent mean field approach, with 74. independent of k (as required by the assumption No. 
(3) from our list), we must have <fio = 0. To set 4>q — we will now use our freedom to adjust thepparameter Rq. As 
discussed in the Section II, we are free to vary this parameter, since it does not affect the dynamicaHJ. Thus we arrive 
at the following condition on Rq: 

0o = ^Co 2 J log(X k jR )K Q (X Jtk /Co)d 2 r 3 = 0. (51) 

After some simple but lengthy transformations, given in Appendix C, Eq. ( pl| ) becomes the following approximate 
expression for log Rq : 



1 \ / 1 



(52) 



So \log(Ro/R)/ \k (R/( )/ y ' 

Note that once we find Co, we can use this equation also to determine Rq. This will not be necessary since Rq drops 



log(R /R)/ \K (R/( 
Substituting Eq. (|52|) in Eq.ji^) yields a closed equation for Co, with no dependence on Rq 

1 \ _ / 1 



out of the dynamic equations, as shown below. Equation ( J53| ) is almost identical to Marqusee's expression for the 
screening length. In principle we can now proceed as planned; for any droplet configuration evaluate Co, calculate the 
mean field approximation to the matrix 0, 

4> itj w K (X itj /Co) fc , fc w (54) 

substitute in Eqs. (38) and d3^) to get L^ 1 and use it in the dynamic equation (p|). There is a problem with doing 



this, though. An important property of the exact L 1 is that it satisfies the sum rule ( |27| ) and hence the total area 
of all droplets is conserved by the dynamics (see Eq.(p9|)). Since Eq.(|54|) is an approximation, we have to check the 
extent to which the sum rules are satisfied by our approximate <fi. Substituting (|54|) in ( |36| ) yields 

f-i _ 1 f-i = K (X i:j /Co) 

i ' i RflogiRi/RoY *J RtlogiRi/Ro) RjlogiRj/RoY 
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Note that L i j are small at distances JQj 3> Co! hence the parameter Co should be interpreted as a screening length. 
At short distances the solution of Eq.([47|) behaves like the lowest order approximation (i.e. stopping at n — 1 in the 
expansion (p4)). Using these expressions in the sum rule (|27|) gives 



R t log(Ri/Ro) 



1+ E 535755^/*') 



(56) 



When we approximate the sum by an integral, in the spirit of our mean-field approach, it becomes 

£ M^)^^'/^) ^ -^o a /*o(r/C)d-r = -1 

and the right hand side of ( p6| ) vanishes, so that in the mean-field limit (i.e. for vanishing area fraction; see Appendix 
D) when our approximations become exact, the sum rule is indeed satisfied. If, however, we use Eq.(|55|) for practical 
calculations at a non vanishing area fraction, the sum rule (and hence area conservation) will not be satisfied rigorously. 



This problem can be overcome rather simply in the following way. Adopt the mean-field approximation (55) for 
the off-diagonal elements of and use the sum rule ( p7| ) to determine its diagonal elements: 

The result for both diagonal and off-diagonal elements of L _1 can be summarized as 

, -1 y K (X^/(o) .,. . 

£-1 _ J RJlog(H ( /flo) ^J(^) logCKj/flo) J / 57 n 

■ KojXij/Co) otherwise 

RtRj logCfli/flo) log(fl 3 /fl. ) omeiwise. 

With this approximation for the sum-rule is of course exactly satisfied in every step and area is conserved. In the 
mean- field limit the sum in Eq.(^) yields —1, and the expression for L~\ obviously reduces to the naive one, given 
byEq.©. 

Using Eq.(]57|) in Eq.(^8|) we obtain now a rather elegant expression for Ri, the rates of change of the droplets' radii: 

r -Vr 1 - 1 v goC*W<b) /j i_\ f ^ 
1 fa" RilogiRt/Ro) ^ togfa/Ro) \Rj Ri)- [> 

For a given droplet configuration we have now an explicit expression for the dynamics of the droplets' radii (note that 
the values of the parameters Rq and Co are also determined by the configuration). As a final "cosmetic" adjustment, 
we use Eq.(^2|) to replace logarithms by Kq (thereby eliminating Rq from the dynamics), yielding the equation we 
used in our numerical study: 

r = 1 V K °( Xi M (1_ 1_\ r591 

1 RiKoiRi/Co) fa, KoiRj/Co) \Rj Rj ' V ' 

To complete the treatment we now justify our naive approach, discuss the regime in which our approximation is 
expected to hold, and check whether various assumptions that were made are self-consistent. First of all we assumed 
that Rq can be chosen so that Co > 0. Now we can see from Eq.(|5^) that the "optimal" Rq ~ Co, and since in the 
mean- field limit Co 3> Ri this means that Rq ^> Ri for all i, so that by Eq.([46|) Co is indeed positive and our approach 
is self-consistent. 

The central point of our approach was replacement of the sum Eq.(^TJ) by an integral: 

g ^JRq ^/Ro)^ k - n (i^/^) / log^/Wfe)^- (60) 

3T^i>k 

It is intuitively clear that such a substitution is valid as long as N^, the number of droplets in the screening zone is large. 
Now, after having obtained the solution 4>, we can indeed show that the condition 3> 1 holds when log(<y9 _1 ) 3> 2n 
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(see Appendix D). Therefore the conditions for validity of the mean-field result lead to a self-consistent theory only 
in the limit of very low area fraction, tp — > 0; even for ip ~ 0.001, there are only a few droplets in the screening zone. 

For the experimental value ip — 0.13, Co is about the mean droplet radius, which makes our approximation completely 
invalid. Such a small value of £o appears because above we have neglected all correlations between the positions of 
the droplets. In particular, it turns out that each of the droplets is surrounded by a depletion zone, from which all 
possible neighbors are excluded (at least, the distance between the centers of any two droplets must exceed the sum 
of their radii). In practice (as seen from our numerical solutions - see below) the mean diameter of the depletion 
zones d is approximately 2.2(R). The corrections provided by including the effects of the depletion zones are discussed 
in Appendix B. As is shown there, taking the depletion zones into account does not change the previously obtained 
expression for (|l6|), but the mean-field (q is replaced by a larger screening length (£ sc = 2.73(R) for ip = 0.13), which 
becomes comparable to the nearest neighbors' distance. Although this improvement is not sufficient to justify our 
mean-field approach, the numerical results (obtained with the effect of depletion zones taken into account) presented 
in the next section are in a rather good agreement with the experiments that were performed at ip = 0.13. 



IV. THE ALGORITHM AND RESULTS OF THE SIMULATIONS. 



In this section we describe our numerical procedure and the results of our simulations. Since we would like to reach 
the scaling state with a sufficiently large number of droplets (to reduce the effect of fluctuations), we must start from 
initial states with very many droplets and allowing this large system to evolve for long times. In order to do this 
with reasonable computational resources one has to resort to approximations that accelerate the numerical procedure. 
Rewrite equations (^9|) in the following form: 

^ = Va ^j2 v ^ ( 61 ) 





k a kp \Rp R a J 

where we used the abbreviated notation 

K a%p = K (X a ^/Csc) k a = K {R a /U) S a = R 2 J2. 

The screening length that appears in the argument of the MacDonald functions is ( sc — 2.73R, as follows from the 
analysis of the depletion zones presented in the Appendix B (see Table I at ip = 0.13). 

Say we wish to integrate (j6l|)-(|62|) naively. For each time step we have to evaluate all N(t) velocities V a , with each 
velocity given by a sum of N terms, i.e. N op « N 2 operations per time step. If we have initially N droplets, and 
wish to reach the late stages with a few % surviving, we would have to perform N steps ~ N time steps of integration, 
as was done by Yao et allla. This is so because if the time step is greater than the life time of some droplet, its area 
will become negative at the end of the time step, which leads to an increase of the area fraction of the surviving 
droplets. Therefore, apparently in order to guarantee exact area conservation one must eliminate each shrinking 
droplet separately, restricting the time step by the next vanishing. At the same time, the detailed evolution of the 
smallest droplet is absolutely unimportant for us. Rather we would like to choose r to be not smaller than it is 
needed to ensure that the relatbpichange of (R), the mean radius of the droplets, is small during one time step. For 
this reason, Akaiwa and MeirontZI simply removed the droplets with R < e(R) with e = 0.1 before performing the 
time step. Let f s be the fraction of the droplets that are eliminated in each time step (it is kept constant, to a good 
approximation, by fixing e). Then the number of steps needed for a run is reduced to N e t eps ~ (1/ fs) log N. However, 
their procedure still did not guarantee precise area fraction conservation^-, when the number of droplets got reduced 
by a factor of 5 (from 100 000 to 20 000), the change of tp was about 2%E1 

We require our algorithm to conserve area fraction exactly and, nevertheless, to reduce both N steps and N op . We 
achieved N op ~ N and N steps ~ (l// s ) log A; that is, our scheme requires 0(N log N) operations for the entire 
evolution. 

To reduce the iV-dependence of N op we note that a droplet j3, whose separation from a exceeds considerably the 
screening length, will have only a small contribution to V a . As will be shown below, only droplets from the first layer 
near a have a significant contribution and therefore only the neighbors of a are included in the sum V a = v a ,(i- 
Therefore we have N op ~ N (replacing N 2 ). 
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To reduce N s t eps to N ste ps ~ (l//s) logiV, conserving the area fraction, we should treat accurately the vanishing of 
the small droplets, instead of simply removing them before performing the time step, as is done in ref£3. Putting this 
into practice requires care, however. Let us consider the system of droplets at two consecutive times to and t\ — t + r. 
We call all droplets that survive to t\ "large" and denote them by indices i and j, whereas droplets that vanished 
before t\ are called "small" and marked by indices n,m. With this notation Eqs. (|6l|)-(|62]) become 

3 rn 



Then a small droplet will live for time 



n j 

We choose r to be small enough so that the density of small droplets is small; hence the typical distance between two 
small droplets exceeds considerably the screening length. Therefore to an excellent approximation two small droplets 
do not interact (i.e. K mjTl w 0) and the first term in ( |64| ) can be neglected. Thus Eq.(|64|) becomes 

= ^Z V m.J =V m . (65) 

j 

t m = -S m /V m . (66) 

Turning now to integrate Eq.(|6^) we note that the term Wj )TO , which represents the contribution of the fast droplet 
m to the growth of the slow droplet i is actually present only during the interval t m and not the entire r. Hence we 
must use 

Si(t + t) = Si(t) + ^ v i,3 T + E v i,mtm 
\ 3 J m 

which can be rewritten, using Eq. (|66|) , as 

Si ( t + r) = Si(t) + (»... + \ E S ^r I T - (67) 

\ 3 m m J 

This is, in fact, the discrete-time form of the differential equation 

3 m 

This equation has a very simple interpretation: the i-th large droplet obtains from each of the small droplets a part 
of the latter's area, proportional to the strength of the interaction between the droplets. Eq.(|67]) is the main working 
formula of our algorithm. The two terms in the parentheses correspond to (i) redistribution of the material between 
the large droplets and (ii) absorption of the material that leaves the small droplets onto the large ones. Using Eq. (|67|) 
eliminates the fast scale of the dynamics; the time scale r is to be chosen so that the assumptions made above are 
indeed satisfied. Our scheme is put to practice by first choosing a convenient fraction of small droplets, which is then 
kept (approximately) fixed throughout the run (typically we used f s ~ 0.003, that corresponds to removing ~ 80 
droplets at each time step when the system contains 28000 droplets). 

Finally, we should take into account the shift of the droplets. For that one can use either the exact relation 

dt 2 2.| X ..|| X .. r ^ 

where the charges qi — RiRi are evaluated at each time step, or the heuristic formula 

dn _ Kjj f 1 1 \ x itj 

dt j^kiki \Ri Rj \Xi,j\* ( ' 

discussed in the Introduction and Appendix A. Then we proceed according the following steps: 
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1. Pick R s such that the number of droplets with R n < R s is approximately N s = f s N. (In practice R S /(R) is 
kept fixed - in the scaling state this is the same as fixing / s ). Identify these as small droplets and calculate all 
V m using Eq.©. 

2. An estimate r^ -* for the time step is determined, as the interval during which all small droplets vanish; 

r<°) = max{i m }, 

m 

where t m is determined by Eq. (|66|) . 

3. Calculate the velocities of the large droplets, Vi = Y^j v i,j (summing over neighbors of i only - hence this step 
takes ~ N operations). Note that some droplet i which has been classified as "large" may, in fact, disappear 
during the interval t' ); this happens if the velocity is such that Si < — ViT^°\ All such droplets, if any, are 
collected and reclassified as "small", a new value for r is determined and the velocities (both V m and V.) are 
recalculated. Usually one such iteration suffices to reach a self-consistent classification. 

4. Integrate Eq.(|6S|) one step r; that is, calculate new areas according to Eq. (|67j) . 

5. Calculate the shifts of the droplets according to either Eq.(|S9|) or Eq.(f70|) and repeat the procedure. 

Since the relaxation time needed to reach the scaling state depends on how close the initial configuration is to this 
state, one would like to choose the initial configuration reasonably close to it. We prepared our initial state as follows. 
We generated a set of N = 50000 droplets with a distribution of radii close to the expected one (determined from 
short preparatory runs). We scattered them over the plane at random but so that each droplet is surrounded by 
a small depiction zone (free of other droplets), to mimic the effect of correlations that appear in the real evolving 
system. For this initial configuration we started a rather large time step r (and run the dynamics without moving 
droplets) yielding only an approximation to the true dynamics. This preliminary relaxation went on till the number 
of droplets was reduced to N = 28000. At this point a smaller time step was selected, the droplets were allowed 
to move and we started to take measurements. All the results presented below are averaged over 8 runs in order to 
reduce statistical fluctuations. 

In our first simulations we did not take into account the shift of droplets. As is shown in Appendix A, this frozen 
droplets approximation might be good even for fractions as large as tp = 0.13. In order to test this, we executed 
runs with frozen droplets and measured the fraction f c of droplets crossing each other, as a function of the number 
of droplets in the system. This data is presented in Fig. 1. It shows a considerable growth of f c with time (from zero 
at N = 28000 to 10% at N — 1000) and there is no tendency towards stabilization. This proves that the droplets' 
motion has a considerable effect on the dynamics of the system and, hence, the model of frozen droplets is invalid for 
such large ip. 

Taking the droplets' shift into account improves drastically this situation. We used both Eq.(|6^) and Eq. (|70|) (called 
models A and B respectively; see below) for the droplets' motion. Fig.l shows that both models A and B give rise 
to much smaller values of f c and that it saturates at long times. Interestingly, while there is no essential difference 
between these two models, the heuristic model B exhibits slightly lower values of f c than model A. In order to test the 
convergence of the sum in Eq. (|69|) we executed runs with different numbers of terms in the sum taken into account; 
namely, we varied the size of the summation box from b — 4.1^ sc to b = 5.8£ sc and found no noticeable difference in 
the behavior of /. Thus we conclude that we have achieved proper convergence of Eq.(|69[). 

We now present results of our simulations, performed at the same relative area fraction, ip — 0.13, that was studied 
experimentally!]. We measured the position correlation function, G(r); that is, the mean number of the droplets whose 
centers lie within an annulus [r, r + dr] around the center of a given one. Fig 2. presents G{r) obtained running model 
B at three consequent moments of time, when the system contained 3000, 2000 and 1000 droplets respectively (each of 
the curves is averaged over 8 runs). We see that the three curves coincide up to their fluctuations, which proves that 
the scaling state has indeed been reached. The corresponding experimental data are also presented in this figure; a 
good agreement is seen. Fig. 3 compares the position correlation functions obtained by the models A and B (averaged 
over time in the scaling state). This proves that our heuristic formula works perfectly. The experimental G(r) has a 
noticeable maximum at r w 4.7, while for our curves the maximum is smaller; at the same time our curve is rather 
close to the result of Masbaum (see Fig. la irurefo), that has a small peak as well, which can not be distinguished 
clearly from the fluctuations. Note that in ref£-l the data were taken with TV ~ 800 droplets, whereas our data is well 
averaged (effectively it corresponds to a system of about 50000 droplets in the scaling state), and the small peak is 
definitely observed. 

Analysis of the distribution of the droplets' radii at three consequent moments of time, where the system contained 
3000, 2000 and 1000 droplets respectively, also shows that the distribution achieved its stationary form. No clear 
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difference between the models A and B was seen. The droplets' radii distribution in the scaling state, averaged over 
time, as obtained from our simulations (model B), is shown in Fig. 4 together with the experimental results. A certain 
discrepancy is seen which, however, does not exceed by much the statistical errors of the experimental points. 

Finally, we measured the charge correlation functional, that contain more detailed information about the system, 
defined as follows. For a charge qi calculate Q + (r), the total amount of similar charge as qi within an annulus [r, r + dr] 
around fi and define the function 

9+{r) = (qiQ+{r))- 

Similarly we define g~(r) in terms of the opposite charges. These two functions, as obtained by simulation of the two 
models, are presented in Fig. 5 together with the corresponding experimental data by Krichevsky and Stavans. The 
agreement can be characterized as excellent. Again, there is no noticeable difference between the results of models A 
and B. 

One-should notice that our definition of the charge correlation functions g±(r), is not precisely the same as that 
of ref.ll Krichevsky and Stavans smeared a droplet's charge on its perimeter before calculating Q±(r), whereas 
we assigned a droplet's charge to it's center. We believe, however, that the smearing of the charge on the droplet's 
perimeter is no more than a way of smoothing the data and there is no essential difference between the two definitions. 

V. SUMMARY 

Ostwald ripening is the coarsening process, observed during the late stage of the evolution of a two-phase system, 
where the droplets of the minority phase exchange material by means of diffusion. This process leads towards a 
scaling state in which the characteristic length scale grows with time according to the scaling law R ~ t 1 / 3 , while all 
the statistical properties (such as the droplets' size distribution, position correlation functions etc.), once rescaled, 
remains fixed. 

The problem of calculating these characteristics has been studied in a number of detailed numerical simulations, 
that take into account all the complicated interactions between the droplets, mediated by the diffusion field. These 
calculations, although being exact, do not contribute much to a qualitative understanding of the importance of 
different components of the interaction between the droplets. On the other hand, analytical mean field treatments 
neglect all spatial effects and seem to be oversimplified. 

The aim of this paper was to construct and test a "minimal extension" of the mean field approach, that will take 
into account spatial effects, keeping only the simplest interaction between the droplets. We calculated analytically an 
approximate form of these interactions using a mean-field approach. Only pairwise interactions between the droplets 
were preserved. We proposed a very efficient numerical algorithm, which allows us to follow the evolution of tens 
of thousands of droplets. We tested our approach by comparing its results with the experimental data and found 
surprisingly good agreement at a relatively large value of the minority phase area fraction, ip = 0.13, where our 
approach to the interaction between the droplets was not expected to work. 

Trying to find the simplest model which reproduces the experimental data, we examined the importance of a number 
of effects. Our findings are summarized as follows: 

1. Depletion zones have a considerable effect on the screening length, increasing it from £ = 1.88/2 to £ = 2.73R. 
At the same time, as shown in Table I, the presence of the depletion zones almost does not affect the functional 
form of the pairwise interaction between the droplets. 

2. Our approach is based on the assumption of circular droplets and monopole + dipolc approximation for the 
diffusion field. The obtained agreement with the experiment indicates that at the values of ip studied higher 
multipoles can be neglected. 

3. Even though inclusion of the depletion zones increases the screening length £ too little to provide formal validity 
to our approximation, our simulations show that it works even for ip = 0.13. 

4. The effect of the droplets' motion is very important for large area fractions, although, formally, it could be 
regarded as adiabatically small compared to the droplets growth. 

5. The expression that determines the shift of the droplets requires summation over a large number of droplets. 
We propose a much simpler heuristic formula (containing a sum over the nearest neighbors only), that gives 
even better results than the exact one. 
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An advantage of our method is its computational efficiency; we are able to choose time steps no smaller than required 
by physical reasonability, eliminating a large number of droplets at each step. At the same time, the total area of the 
droplets is conserved exactly at each time step. This makes our approach useful for extensive studies of the Ostwald 
problem. 
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APPENDIX: A ANALYSIS OF THE FORMULA FOR THE DROPLETS' SHIFT. 



The droplets' shifts are determined by Kq.(pO|): 

~dt 



Xi.4 \Xa 



(Al) 



First of all, using this formula we can show, that in the low area fraction limit the shift of the positions can be 
neglected. Indeed, T s h, the characteristic time of a shift of a droplet's center of mass is given by 



1 



dX, 



dt 



The characteristic growth time of the droplet is (see Eq.(12)) 



gr 



1 



dRi 



dt 



RT 



Comparing Eqs.(A2) and 



we see that 



-l R 2 i 



.f 



9-r X 2 - 
1,3 



9r 



(A2) 



(A3) 



(A4) 



That is, for small area fractions the motion of the droplets' centers is adiabatically slower than their growth. Conse- 
quently, one can neglect the droplets' motion and the system is characterized only by the dynamics of the droplets' 
radii as determined by Eq.([D|)- (p^). Formally, one can expect this approximation to be valid even for ip = 0.13 
(used in the experiments by Stavans and Krichevsky) and we have tried it in our work (see Section IV and Fig. 1). 
Our simulations have shown, however, that neglecting the droplets' motion gives wrong results and, therefore, the 
dynamics of the fi has been taken into account. 

Secondly, note that a rough estimate of the sum on the r.h. of Eq.(Al) indicates that it exhibits bad convergence 
properties. Assuming qj to be uncorrelated random variables with zero mean we get 



dr.; 
dt 



y— 



(g 2 )iogi s , 



(A5) 



where L s is the size of the system. A more accurate treatment implies, however, that the sum does converge. According 
to Eqs.®-® 



E K 



where 



1 



1 



1 



k m hj \ R m , Rj 



and we use the abbreviations Kj^ = Ko(Xj^/C sc ), kj = Ko(Rj /( sc )- Then Eq.(Al) becomes: 

dXi \ ^ (j) Xij \ ^ \ ^ ri( m ) *>J 



(1/2)- 



dt 



\X 



h3 I 



E E Q) 

j (^i) m( =£j,i) 



\x^r 



(A6) 



where 



Qf=h \ 



denotes the part of the charge on the i-th droplet that is induced by the j-th one. The first term of Eq. ( |A6[ ) determines 
the shift of the droplet due to the direct material transfered between this droplet and its j-th neighbors. This sum 



16 



converges very well because decreases exponentially with distance X^j. The second term accounts for the effect 
of the redistribution of the material between the j-th and m-th droplets. Since the internal sum (on m) contains a 
short-range factor Kj t7n , it is actually over a(x (-box around the j-th droplet. In this double sum, each term 

Q{m) XjJ 



can be paired with: 



X, 



i . in 



2 • 



Since = — Qm , these two contributions can be considered as a dipole. Thus, each ( x (-box represents a dipole 
P (randoml y di rected) with the dipole moment P ~ q(, where q is the characteristic scale of Qm ■ Thus, the second 



term in Eq.(A6) is now estimated as: 



second term ~ £ - 2{P ~ ]T -f-^. (A7) 

The mean of this expression is zero while the mean square deviation is given by 



which converges to a finite value. Thus, when calculating the sum in Eq.(Al), we can restrict ourselves to only several 
nearest layers of neighbors. 

Finally, one can use an even simpler heuristic formula for calculating the shift of the droplets. The meaning of 
Eq. (|A6|) is that the motion of z-th droplet has two sources. The first is the material transfered between this droplet 



and its j-th neighbors (the first term of Eq.(A6)). The second is due to redistribution of the material between the 
surrounding j-th droplets them selves (the second term). Although these contributions arc of the same order, in 
our case we have a reason to drop the second term (although it does not simplify computations, it does makes the 
model physically simpler). The shift of the droplets has a noticeable effect only at relatively large fractions, where 
the interaction between the next nearest neighbors is considerably suppressed. Then, for a fixed configuration of the 
nearest neighbors of the z-th droplet w e ca n vary the configuration of its next nearest neighbors. This manipulation 
will not affect the first term of the Eq. (|A6| ), while it will reduce its second term. Thus, in the mean field spirit of our 
model, we can average the shift velocity of the z-th droplets over various configurations of its next nearest neighbors. 
Thus, we finally get: 

<l*; 2 V^fl-i)A ? . ( A9) 



dt k-^kj V Pi'i Pvj J I Xl.^ j 1^ 

Although the approximation leading to this for mula is not based on a rigorous expansion, our simulations show 



that it works as well as the more rigorous Eq.(Al). At the same time it is much more economic because it requires 
the summation only over the nearest neighbors. 
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APPENDIX: B DEPLETION ZONES 



The simplest possible improvement over the mean-field approximation can be obtained by including some of the 
correlations between the positions of the droplets' centers. In this Appendix we will take into account the simplest 
manifestation of these correlations — the so called depletion zones. These are the regions around each of the droplets, 
from which all possible neighbors are excluded (by geometrical steric constraints - the distance between the centers 
of any two droplets must exceed at least the sum of their radii) . 

The exact equation for the diagonal elements of the matrix cf> is given by Eq.(^l]), whereas for the off-diagonal 
elements it is Eq.(39) 



1 



. log Rj/Bi 



■]og(Xi t j/Ro)(j)j t k + 4>i,k = -j k logX^ k /R . 



(Bl) 



This has been replaced in our mean field approximation by the following integral equation for the smooth function 
cf>(x) (see Eq.©): 



1 

2^' 



Co" 2 / log(X IJ / J R o )0(^, fe )d 2 r, + <f>(X ilk ) = -jlogX iik /Ro. 



(B2) 



Clearly, by using the integral of Eq.(B2) wc implicitly assume that the distribution of fj, the positions of the centers 



of droplets j, are independent of their distance from the one at f. This homogeneity assumption may serve as a 
reasonable approximation as long as X, the mean distance between neighbors, is much greater than d, the typical 
radius of the depletion zones. When, however, the density increases to the extent that X ~ d, the inhomogeneity of 
the distribution of the j-droplets around the fixed droplets i and k can no longer be ignored. One should emphasize 
that other effects, such as correlations between different "charges" and between the droplets' sizes may also be of 
importance and can also possibly affect the elements of the inverse matrix we are calculating. Nevertheless, here we 
take into account only the correlations that ensure that the areas of two neighbors do not overlap. In our actual 
numerical calculations even this is done only approximately, by introducing a uniform sized depletion zone, neglecting 
its fluctuations as well as dependence on the droplets' radii. 

Correlation b etw een the positions of the droplets can be incorporated in our approximate treatment by replacing 
the sum in Eq.(Bl) by the following integral operator: 



Sit 



Z7T 



\og(X h] / R Q )P(n,r J ,r k ) ( f,(X 3 .k)d 2 r 1 



(B3) 



where P is the probability of finding a droplet centered at the point f = fj, given that there are droplets at f = f 
and f = f k . We approximate this probability (three-point correlation function) by representing it as the product of 
pair correlation functions: 



P(n,fj,f k ) = g(X itj )g(X jtk ), 



(B4) 



where g(Xij) is the probability of finding a droplet j at distance Xij from the center of droplet i; this function is 
normalized such that g{r) — * 1 at r — ► oo. Within this approximation Eq.(Bl) is replaced by 



Sint4> + 4>{ X i,k) 



5 k, k 



\og{R k /Ro) 



\ogXi^/R 



(B5) 



where now we have 



W = -^Co / log(X lJ /R a )g(X lJ )g(X J ^(X J , k )d z r 1 . 



(B6) 



In order to solve equation (B5) we first try to bring it as close to the form of Eq.(B2) as we can, and then use the 
same method of solution as was used there. To this end we first rewrite the expression for SVnt as follows: 



S int <t>= -^Co~" / log(X itj /R )g(X itj )c} > (X jtk )cC 2 r j 



2tt 



Co" 2 / log(X id /R )g(X h] )[l - giX^MXj^cPrj 



(B7) 
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The pair correlation function g(X) vanishes for short distances, X < d, and g « 1 for large X. Therefore, we get 
non-vanishing contributions to the second term in Eq.(B7) only when Xj k < d. For Xi )k 3> d, we have in this region 
g(X itj ) w 1 and log(X itj / Rq) w \og(X^ k /R ), so that 



5 int 0«-— Co" 2 / log(X i , J / J R ) 9 (X i , J )^(X j>fe )d 2 rj 

Z7T 



2tt 



Co"' \og(X lik /R ) / [1 - ./( .\V / .ro(.\V / .) ( / J r ; . 



(B8) 



The last formal step we take is to express the first integral here as the sum of two terms, using —g = [1 — g] — 1, 
which leads to our final expression for Si nt 4>: 

1 



+ ^Co" 2 / log(X iJ /i? ) [1 - <?(^)] 0(^,*)«?r j 
1 

~ 2tt' 



C -^log(X i>fe /iio) / [1 - g(X j>k )]<f>{X jti )d l r j 



(B9) 



Our basic equation ( |B5| ) takes now the form 
1 



2tt 



—Co" 2 / log(X^/i? )^, fc )^r i + ^(X i , fe ) + — Co"" / F(XijWX jJt )#r 



where 



and 



-7fc log(X ijfc /i? ), 



7fc = 1 + 



+ ^Co"" l(l-g(X J , h ))<l>(X j , k )d'r J 



log{R /R k ) 2tt 



F(X)=log(X/i? )[l-9(X)]. 



(BIO) 



(Bll) 



(B12) 



Note that equation ( BIO ) is very close in form to Eq.([B2|); the only significant difference being the appearance of 
a new term - the integr al ov er the function F, which contains explicit ly th e correlation function g. Clearly when 
g(X) e 1 we recover Eq.(B2). As was the case there, the solution of Eq.(BlC) is symmetric if and only if ^ k does not 
depend on k, 7^ ~ 7, which is achi eved by imposing cf> k-k = 0. This leads again to a condition that determines the 
value of the parameter Rq (see Eq.(B21) below). Although now 7 7^ 1, (see Eq.(Bll)), its value is unimportant for 
the simulations since it can be absorbed in the definition of the time scale in the dynamic equation. Hence below we 
will set 7 = 1. 

We solve Eq.( B10] ) by first applying the operator to the equation, which yields 

(B13) 



- Co" y(X ilk ) + V^(X a ) + ^Co""V^ / F {X it j )<fi(Xj t fe ) d tj = -2^S(X iik ). 
Denote the Fourier transform of F(X) by 

= Jo rdrl ° g (j^) ^~ 9i r )]M<l r )- 



(B14) 



Note that the integral converges, since g(r) — > fast for large r. Upon Fourier transforming Eq. ( B13 ) becomes 

- C - 2 <K<7) - q 2 kq) - (o 2 q 2 F(q)4>(q) = (Bi5) 

yielding the solution for <fi(q) (from now on we set 7 = 1) 
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2?r 1 



(B16) 



Once 4>{q) has been evaluated, the function <fr(X) is obtained by the inverse Fourier transform. ^From this point on 
we work with the specific (step-function) form of g(X): 



X < d 
X >d 



(B17) 



For this form of g the integral in Eq. (B14) can be calculated analytically, integrating by parts and using the identities 
[xJi(x)]' — xJq(x) and [Jq(x)]' = Ji{x), to get 



F{q) = q-^-qdJx^logiRo/d) + [1 - J (qd)}}. 



(B18) 



We now substitute this in ( B16| ) and perform the angular integration in the inverse Fourier transform, to arrive at 
the following expression for <fi(r): 



0(r) 



■Mq r ) 



o q 2 + Co - [Co log(i?o/d)]gdJi (qd) + Co' [1 - J (qd)] 



qdq. 



(B19) 



We cannot perform this integral analytically, but it is simple to evaluate by numerical integration. Once this has 
been done, the parameter Rq can be found by imposing the condition (f>k k — 0, which takes the form (analogous to 
Eq.©) 



Ok,k 



J_ C - 



log ( — ) g{r)(j){r)d 2 r 



which becomes the following equation for Rq: 

log Rq = 



XT log(r)(p(r)rdr 



IT ${r)rdr 

Once Rq has been determined we can evaluate the parameter Co, using its definition 

1 



C z = 27rn< 



log(iV^)' 



(B20) 



(B21) 



(B22) 



For a given droplet configuration one can now proceed to find (j>(r) by taking the following steps. First, determine 
the radius of the de pletio n zon e d from the measured correlation function of the droplets (such as Fig. 3). Then 
solve Eqs.( Bf 9 , B2f ) and ( B22 ). This can be done iteratively by initializing the proced ure u sing the results of the 
correlation- free t heory for <j>,Co an d ^?o- The next iterate of cf>(r) is evalu ated using Eq.(B19); this new <fi(r) is then 
used in Eq.( B21 ) to yield the new value of Rq; this, in turn, is used in Eq.( B22 ), together with the distribution of the 
droplets' radii, to yield the new iterate of Co- The procedure is then repeated until convergence is reached. In practice 
we used a few simplifying steps, which made computation faster without significantly altering the result. The first 
simplification consists of setting 



d = k(R) 



(B23) 



for an entire run, instead of determining it from the correlation function at each time step. We found that values in 
the range 2.0 < n < 2.8 fit the correlation function quite well. The numerical results were obtained using the fixed 
value k = 2.15. A second time-saving simplification is to use the approximation 



Co 



2im 



1 



log(i? /(i?)) 



(B24) 



to calculate Co, instead of evaluating the average (l/log(i?o/-Ri)) over the measured distribution of droplet sizes. 
Recall also that the density of droplets, n, is related to their relati ve ar ea fraction, ip, by n — ip/Tr(R 2 ). In the spirit 
of the previous approximation we replace (R 2 ) « (R) 2 , so that Eq.(B24) becomes 



c - 



2^ 



(R) 2 log(R /(R)) 



(B25) 
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Using Eqs.(B23) and (B25) simplifies our numerical scheme significantly. To demonstrate the effect of including the 
depletion zones in the calculation we present in Table I the results of calculations performed in the scaling regime at 
four different values of (p. The various quantities presented were obtained as follows: Rq and Co are (simultaneous) 
solutions of Eqs. (p!9|) , ( p2l| ) (with d = 2.15(i?) in both) and ( |B25| ). (™ f 
obtained by setting d = 0. The screening length C 



is the value of the screening length as 
is the first moment of the function </>(r), 



c 



2 

scr 



4>{r)rdr 



(B26) 



calculated by numerical integration of Eq.(BlE). 

As we see, for very small fractions the screening length is less than Co • This is surprising, since we expected 
inclusion of the depletion zones to increase the screening regions. For the larger fractions, however, we indeed see 
that the effect of including the depletion zone changes sign: the renormalized screening length becomes greater than 
the mean field result, as was expected. 

Another important observation we should make is that as tp increases, C decreases and becomes comparable to the 
average radius (in units of which it is given in Table I). At such </?, C obviously cannot be interpreted as a "screening 
length" , because C « (R) means that there are no droplets in the " screening zone" . 

Another question that we studied in detail concerns the extent to which introducing the depletion zone affects 
the function 4>(r). Clearly, when we are not in the limit of very small ip it is no longer given by the MacDonald 
function. Table II contains </>(r), as obtained at <p = 0.13 by the procedure outlined above: setting d — 2.15(R), 
and simultaneously solving Eqs.( B19 B21 ) and (B25). Comparing 4>(r) with the simple mean- field result K (r/C m/ ) 
reveals that the renormalization due to inclusion of the depletion zones leads to significant change of the matrix 
elements 4>(r). On the other hand, comparing 4>{ r ) with the function Ao(r/Co) we see that (f>(r) does not differ much 
from MacDonald's function, provided the renormalized Co is used. Therefore in principle this function can be used 
in calculations as a fairly good approximation for 4>(r). The numerical results presented in Section IV were obtained 
using Co = 2.73, according to Table I. 



APPENDIX: C DERIVATION OF THE CONDITION FOR R . 

The condition (45) on Rq: 



^Co"" / n»KlA' / ,.,./i? ) ,.)/v )( (A' J -. i ./g,)^/ :/ = 



(CI) 



can be simplified by representing 

so that Eq. ( |Cl| ) becomes 

where the constant A has the value 



log(Xij/R ) = log(Xij/Co) +log(Co/i?o) 
B + Alog(Co/i?o) = 



because of the normalization of Kq(x), while the other constant, B, is given by 

B = ^Co 2 J log(A fcj /Co)A (A 3 , fe /Co)d 2 r, 



(C2) 
(C3) 

(C4) 



In order to evaluate B we note that the solution of Eq.(j45|) is given, for any value of Rq, by Eq.([4g|). Therefore we 
can substitute 4>{Xi^) from Eq.(48) in Eq. (f45|) and choose for Rq the special value Rq — Co, to get the identity 



2tt 



Cq 2 / logflf - r'|/Co)A (r7Co)d 2 r' + ^o(r/Co) = - log(r/Co) 



(C5) 



It is trivial to see that in the limit r — > the integral on the l.h.s. of this identity becomes precisely equal to B (as 
given by Eq. (|C4|)). For f-+0we can use the small r limit (E9J) of Kq so that for very small r our identity becomes 
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so that we find 



-.B - log ( — I - C = ],., 



B = log 2 - C 



(C is Euler's constant). Using the values of A and B in Eq.(C2) it becomes 

log2-C + log(Co/i?o) = 



(C6) 



(C7) 



providing a new connection between Rq and Co- This is to be used together with Eq.(f46|) to solve for both Rq and Co- 
Since Eq.(E^) contains log(Ro/(R)), it is convenient to add and subtract log(i?) from Eq.(C7) and to rewrite it as 



log(i? / < R >) = log(2Co/(i?» - C 



(C8) 



If we have (R) < Co the right hand side of Eq.(|C8|) is « K ((R) /Co), so that (|CT|) becomes log(i? / < R >) ~ K (< 
R > /Co)- Finally, we can write, with the same accuracy: 



(C9) 



\og(R /R) 



Ko(R/(o) 



APPENDIX: D THE SMALL PARAMETER OF THE APPROXIMATION. 

In the mean field limit, when i?o ~ Co 3> Ri the quantity 1/ \og(Ro/Ri) does not vary much and one can take it 
out of the sum and replace it by its mean value. Secondly, we see that there is a scale Co such that <f) does not change 
much over distances much less than Co- Let us divide the plane into boxes b of size Co x Co- Then the sum that still 
remains can be rewritten as 



E logprWW** 



(Dl) 



Actually, only the few boxes located near have a significant contribution to the sum; the contribution of all the 
others is exponentially small due to the screening effect. 

The expression being summed does not change much inside each box and therefore, if the mean number of the 
droplets in the box, Nq, is large enough, the internal sum within each box can be replaced by the integral. In this 
case, according to the theorem of large numbers 



^logiXij/Ro^k-n f log(X h] /R )^(X Jtk )d 2 r 3 



1 



This gives rise to the following condition for the validity of our approximation: 

- "Co 2 » i 

But Eq.(^6|) implies that nCo ~ (Stt)" 1 \og(Ra/(R)) and, as we have shown Co ~ Ro, so that 

nC 2 ~ log(Co/<#» 

On the other hand, multiplying Eq.(|4^) by (R) 2 and using the definition of the area fraction ip — Tin{R 



yields log(Co/-Ro) ~ log(v 1 )- This, together with Eqs.(p3-D4) means that 

log(</? _1 ) > 2?r. 



(D2) 

(D3) 

(D4) 
irn{R) 2 
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Ro 


so 


Co 


Cscr 


0.001 


51.35 


45.99 


45.96 


45.922 


0.01 


15.59 


11.77 


12.13 


11.71 


0.05 


8.10 


4.08 


4.73 


4.00 


0.13 


6.14 


1.88 


2.73 


1.94 



TABLE I. Various quantities related to the self consistent determination of <j>(r) (with th e depletion zones taken into 
acco unt) for different area fraction ip. Ro and (o are (simultaneous) solutions of Eqs. (B19,B21) (with d — 2.15(7?) in both) 
and ( B25 ); £ m / is the valu e of t he screening length as obtained in the Section III (neglecting the effect of the depletion zones); 
£ scr is determined by Eq. (B26) with </>(r) obtained by numerical integration. All lengths are given in units of (R). 



r/x 


4>{r) 


Mr/Co) 




.722 


2.884 


2.617 


1.205 


.884 


1.914 


1.773 


.711 


1.045 


1.320 


1.218 


.425 


1.206 


.901 


.845 


.257 


1.420 


.546 


.525 


.133 


1.689 


.297 


.294 


.059 


1.904 


.185 


.186 


.031 


2.012 


.147 


.148 


.022 


2.118 


.118 


.118 


.016 


2.333 


.075 


.076 


.008 


2.548 


.046 


.048 


.004 


2.763 


.027 


.031 


.002 


2.97 


.016 


.020 


.001 


3.031 


.014 


.018 


.001 


3.353 


.0087 


.0095 


.00047 



TABLE II. Interaction function <j>(r) for the case of depletion zones calculated at the area fraction ip = 0.13. It is compared 
with the MacDonald's function with the renormalized screening parameter £o (third column), and with the unrenormalized, 
mean-field screening length (fourth column). Note that the distance r is measured in units of the mean distance between 
neighbor droplets, x = 1/y/n. 
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FIG. 1. Time dependence of f c , the fraction of crossing droplets, as obtained from the frozen droplets model (empty 
squares); from model B, that uses the heuristic Eq. @ for the droplets' shift (full squares) and from model A, that uses the 
exact Eq. d&s|) for the droplets shift (large circles and small circles for small and large summation block sizes: b = 4.1£ and 
b = 5.8£ respectively; see explanations in text). The lines are guides for eye. The growth of f c observed for the frozen droplets 
model indicates that it is not valid for the area fraction tp = 0.13. On the other hand, once the droplets are allowed to move, 
f c decreases to very small values irrespective of whether model B or A is used. 
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r 



FIG. 2. Correlation functions of droplets' positions, G(r), as measured at different moments of time in the scaling state, 
corresponding to N = 3000, N = 2000, N = 1500, N = 1300, N = 1000 droplets present (small ctrcles). The weighted average of 
these runs (the weight of a configuration is proportional to the corresponding number of droplets) is also shown (large circles). 
These results were obtained using model B. Note that all the lines are close to each other, indicating that system has indeed 
arrived at the scaling state. The experimental results of Krichevsky and Stavans, shown for comparison (squares), are also close 
to ours. 
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FIG. 3. Weighted time average (the weight of a configuration is proportional to its number of droplets) of the correlation 
functions of droplets' positions, G(r), (averaged over 8 runs) for models A and B in the scaling state (large circles and small 
circles respectively). The experimental results of Krichevsky and Stavans are also shown (squares). This shows that the two 
models give indistinguishable results that are in good agreement with experiment as well. 
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FIG. 5. The charge correlation functions (see precise definition in the text) for the same (g+(r)) and the opposite (g~(r)) 
charges, as obtained using models A and B, compared with the experimental data of Krichevsky and Stavans (full circles and 
full squares for the same and the opposite charges respectively). The lines are the guides for eye. Our data present averages 
over 8 runs. Note that the fluctuations in this plot are larger than for the position correlation function (see Fig. 2). We believe 
that the difference between the results of the two models is due to the fluctuations. Interestingly, model B seems to be closer 
to the experimental points. 
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The role of the parameter Ro in our theory is analogous to the parameter A in Beenakker's modeled; while being irrelevant 
for the exact problem, it was tuned in order to restore conservation of material, broken in the reduced description. 
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